// $Id: rfits.c 172 2013-06-20 03:34:01Z alan.watson.f $

#include <stdio.h>
#include <string.h>

#include "array.h"
#include "rfits.h"

// oaxaca_rfits() implements a very limited FITS reader. It is only
// intended to read 16-bit integer FITS files without scaling.

array_short_t *
oaxaca_rfits(FILE *fp)
{
  int bitpix;
  double bscale = 1.0;
  double bzero = 0.0;
  int naxis = 0;
  int nx = 0;
  int ny = 0;
  
  char record[81];
  record[80] = 0;
  if (fread(record, 80, 1, fp) != 1)
    return NULL;
  int nrecord = 1;
  while (strncmp(record, "END   ", 6) != 0) {
    sscanf(record, "NAXIS  = %d", &naxis);
    sscanf(record, "NAXIS1 = %d", &nx);
    sscanf(record, "NAXIS2 = %d", &ny);
    sscanf(record, "BITPIX = %d", &bitpix);
    sscanf(record, "BSCALE = %lf", &bscale);
    sscanf(record, "BZERO  = %lf", &bzero);
    if (fread(record, 80, 1, fp) != 1)
      return NULL;
    ++nrecord;
  }
  
  while (nrecord % 36 != 0) {
    if (fread(record, 80, 1, fp) != 1)
      return NULL;
    ++nrecord;
  }
  
  if (bitpix != 16 || 
      naxis != 2 || nx == 0 || ny == 0 || 
      bscale != 1.0 || bzero != 0.0)
    return NULL;
  
  array_short_t *a = alloc_array_short(ny, nx);
  
  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix) {
      unsigned char b[2];
      if (fread(b, 2, 1, fp) != 1) {
        free_array_short(a);
        return NULL;
      }      
      unsigned short u = (b[0] << 8) + b[1];
      short i;
      if (u & 0x8000u)
        i = - (short) (~u & 0x7fffu) - 1;
      else
        i = u;
      a->z[iy][ix] = i;
    }
  
  return a;
}
